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ABSTRACT 

We show that the algorithm proposed by Gauss to compute the secular evolution of 
gravitationally interacting Keplerian rings extends naturally to softened gravitational 
interactions. The resulting tool is ideal for the study of the secular dynamical evolution 
of nearly Keplerian systems such as stellar clusters surrounding black holes in galactic 
nuclei, cometary clouds, or planetesimal discs. We illustrate its accuracy, efficiency and 
versatility on a variety of configurations. In particular, we examine a secularly unstable 
unstable system of counter-rotating disks, and follow the unfolding and saturation of 
the instability into a global, uniformly precessing, lopsided (m = 1) mode. 
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1 INTRODUCTION 

The theory of secular variations of the planetary orbits was introduced by Lagrange over two centuries ago. This approach 
, to celestial mechanics is based on the observation that the orientation of eccentric orbits in a Keplerian potential is fixed: in 
' a coordinate representation, the argument of periapsis co and the longitude of the node Q are stationary. Because of these 
conserved quantities, the time-averaged torque between two planets on fixed orbits is generally non-zero if the orbits are both 
■ eccentric and /or mutually inclined. Thus, if the typical planetary mass is m and the mass of the host star is M*, the inter- 
' planetary gravitational forces lead to both short-term oscillations, having periods of order the orbital period P and fractional 
, amplitudes of order m/Mi,, and long-term or secular oscillations, having periods of order {M-^/m)P, fractional amplitudes 
• in eccentricity and inclinations of order unity, and constant semi-major axis. The large amplitude of the secular oscillations 
means that in most cases they dominate the evolution of the planetary system. 

In the solar system, secular variations are usually computed as power series in the eccentricities e, inclinations i, and mass 
ratios m/Mi, of the planets, which are all sma ll. Current se cular theories employ differential equations of motion containing 
~ 10^ terms up to order (m/M*)^, , and (|Laskaij|l986l '). 

Expansions of the equations of motion in powers of eccentricity and inclination work well for the orbits of planets in 
the solar system (median eccentricity and inclination 0.05 and 1.8°), but not so well for orbits with larger eccentricities 
and inclinations, such as those of asteroids (median eccentricity and inclination 0.14 and 7°), extrasolar planets (median 
eccentricity 0.2), the stars that surround the cent ral black holes found in many galaxies, and planets or binary stars executing 
Kozai oscillations excited by a co mpanion star (jlnnanen et al.lll997l : iMazeh et ai1ll997l : iHolman et allll997l : IWu fc Murravl 
l2003l : iFabrvckv fc Tremainell2007l ) . 

The primary goal of this paper is to develop numerical methods for secular dynamics that are valid for large eccentricities 
and inclinations. In so doing we shall restrict ourselves to equations of motion that are first order in the ratio m/Mi,. This 
approximation does not work well in the solar system, mostly because of near-r esonances between the outer planets (e.g., the 
'great inequality' between Jupiter and Saturn; see also iMurrav fc Holmanll 19991 ), but should be more useful in systems where 
the mass ratio m/Mi, is smaller and/or we are seeking a qualitative understanding of the secular evolution rather than an 
accurate ephemeris. 

Although the results we describe in this paper are applicable to any near-Keplerian system, for conciseness we shall 
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refer to the central body as a 'star' and the satelhtes in orbit around it as 'planets', except in ^ where we examine the 
performance of the algorithm with parameters appropriate for stars orbiting a black hole. To first order in the planetary 
masses, the equations of motion for secular dynamics are derived from a governing Hamiltonian that represents the sum of the 
time-averaged interaction potential energies between all pairs of planets in the system. Therefore the central ingredient in any 
numerical method for secular dynamics is the efficient evaluation of the average interaction energy between two fixed Kepler 
orbits, each specified by its semi-major axis a and shape parameters e (eccentricity), u) (argument of periapsis), I (inclination), 
and n (longitude of node). The averaging can be thought of as smearing each planet into a 'ring' in which the mass per unit 
length is inversely proportional to the speed of the planet in its orbit. We call a code that follows the secular gravitational 
dynamics of A'' planets an 'N-ring code' by analogy with N-body codes that follow the full gravitational dynamics. 

The simplest way to calculate the interaction energy between rings labelled by a and (5 is to choose K points equally 
spaced in mean anomaly on each ring and then compute — Gma'mp/{K'^ IS.ij), where rUa and mp are the masses of 

the two bodies and Ai^ is the distance between point i on ring a and point j on ring j3 (e.g.. ICiirkan fc Hopman| [2007') . This 
approach is rather inefficient, both because it requires 0{K^) computations per r ing pair and becaus e it co nverges slowly, 
especially if the orbits have a close approach (minAij ^ maxAij ). For example. iGiirkan fc HopmanI (|2007l ) typically find 
that K ~ 10^'^ is needed for 1 per cent relative accuracy in the potential. 

A better approach, which we pursue in this paper, is to compute the potential from ring a at an arbitrary point r by 
integrating analytically over the mass elements of ring ct, and then to numerically integrate this potential over the locus 
of points r that lie on ring 13. The rather involved analysis needed to do this was accomplished by Gauss and Hill in the 
nineteenth centurjQ. We shall call this approach 'Gauss's method', and the use of Gauss's method in an N-ring code the 
'Gaussian ring algorithm'. 

When he wrote his classic memoir on the force exerted by an elliptic ring on a mass point. [Causd l|l818l ') achieved at least 
three lasting outcomes: he provided dynamics with a vivid illustration of the averaging principle; he introduced the arithmetic- 
geometric mean as an efficient method for computing Legendre's elliptic functions; and he showed how the evaluation of the 
first-order average of the planetary equations can be dramatically accelerated with the help of such functions. iHagiharal (|l97ll ) 
gives a thorough review of Gauss's method, its algorithmic manifestations and its applications. Its main attraction of course 
is the efficient averaging of the equations of motion over the fast, orbital, time scale. The resulting differential equations, in 
the case of the solar system, can be evolved with a timestep of at least 500 yr-over 10'* times larger than the timestep needed 
to follow the unaveraged equations of motion. The main drawback of Gauss's method, as already noted, is that it neglects 
perturbations that are second-order or higher in the planetary masses, some of which are near-resonant terms that, in the 
solar system, can be almost as strong as the secular terms that are first order in t he masses. 

An implementation of Gauss 's method by one of us (J.T.) was used to explore [Kozail (Il962l ) oscillations of the then newly 
discovered planet 16 Cygni Bb (|Holman et al.lll99'<1 ). In most applications, Kozai oscillations arise from the influence of a 
distant companion body. In the case of 16 Cygni Bb, this is a stellar companion to the host star with a projected separation 
of ~ 850 AU, 500 times the planetary semi-major axis. At such large separations the perturbing potential from the stellar 
companion is dominated by its quadrupole moment, and in this case the averaged interaction potential between the companion 
and the planet can be computed analytically. Nevertheless, Gauss's method offers important insight into this problem. The 
reason is that in the quadrupole approximation the averaged interaction potential is axisymmetric in the reference frame 
whose equator is the orbital plane of the companion, even if the companion orbit is eccentric. Thus both the energy and the 
component of angular momentum along the axis of the companion orbit are integrals of motion, and the secular dynamics is 
described by a Hamiltonian with only one degree of freedom. The existence of this second integral is an accidental consequence 
of the properties of the quadrupole potential, and without the quadrupole approximation the integral is not present. Gauss's 
method allows us to explore the two degree-of-freedom problem defined by the full secular dynamics, which for example shows 
that in many cases initially circular planetary orbits are chaotic and allows us to estimate the width of the chaotic zonfl 

What drives much of the current interest in Gauss's method is the wish to study the evolution of a distribution of stars 
within the sphere of influence of the central black hole in a galaxy. In this situation the central object enforces Keplerian 
dynamics, perturbed by the mutual gravitational attraction between the stars. Here then is a natural context for Gauss's 
method. However, in such problems, one immediately faces close encounters between particles, which translate into ring 
crossings. At a crossing, the force between rings chang e s disc ontinuousl;^, a jump that is difficult for numerical integrators to 
follow. Such was the experience of I Ranch fc Tremaind (|l996l ). who used an unsoftened version of Gauss's method to simulate 
resonant relaxation around black holes and found that the speed-up resulting from averaging is virtually cancelled by the 
slow-down caused by ring crossings. 

In N-body simulations, one often resorts to softened gravitational interactions to alleviate the cost of resolving near- 
collisions between particles (and to better approximate coUisionless behavior). It is natural to ask whether Keplerian rings 
that perturb each other with softened gravitational forces can be evolved with the help of a softened variant of Gauss's 

^ It would be better, of course, to integrate over both rings analytically but even Gauss was apparently not able to do this. 

An alternative approach is to expand the potential to octupole order, which also destroys the second integral; see lFord et al. 1 1 I2OOOI) . 
^ Consider a vertical ring with mass per unit length u lying on the z-axis, and a similar horizontal ring lying parallel to the j/-axis and 
passing through the point {x, y, z) = {xo, 0, 0). The force exerted by the vertical ring on the horizontal ring is —2nGa^ sgn (x)x. 
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method. We were pleased to learn that the answer is yes, and we will proceed to tell you how. We start in ^by reviewing 
the underpinnings of Gauss's algorithm, and how they are modified for the softened case. We thereby derive a formulation 
of Gauss's method for softened gravitational interactions. We discuss how to implement the method numerically in Sj3] In ^ 
we describe numerical examples that illustrate its accuracy, its strengths, and its limitations (limitations of secular averaging 
in general), by direct comparison with simulations of the analogous unaveraged system. These tests lead us to !j5l where we 
present simulations of a (linearly) unstable configuration of two counter-rotating annuli, which track the instability to its 
(nonlinear) saturation as a global, uniformly precessing, lopsided (m — 1) mode. We conclude in |6]with a discussion of the 
main features, potential applications, and future improvements of our algorithm. 



2 GAUSS'S METHOD REVISITED 

We argued in the preceding section that to comput e the interac tion potential between two Keplerian rings using K equally 
spaced points per ring required 0{K^) calculations. lGaus3 (|l818l ) showed that the wor k can be red uced to 0{K) calculations 
by integrating analytically over one ring with the help of Legendre's elliptic functions. iHilll l|l882t ) elaborated the arguments 
in Gauss's paper, and described an algorithm for the exact calculation of se cular evolution w hich is accurate to first order 
in the masses and all orders in eccentricity and inclination. Hal phcni (|l888l ). and iGoriachevI (1937) after him, follo wed an 
alternative, geometric, route to elliptic functions, which was later described in tensorial lang uage by iMusenI (l963l ). While 
Gauss's exposition is purely algebraic in fiavor, Halphen exploits the geometry of the cone (with base on the disturbing ring 
and apex at the point where the force is evaluated), to motivate a transformation to the cone's principal axes that makes 
Gauss's insight transparent. While the algebraic manipulations work through effortlessly for softened interactions, it is less 
obvious how to soften things geometrically, since the softening scale breaks a crucial symmetry. We choose to cover the 
algebraic derivation in the body of the text, relegating the geometrical treatment to Appendix El 

The algorithm, softened or not, has two independent components: on one hand, the expression of the averaged equations 
of motion in a fiexible, singularity-free coordinate system, while presuming a means, efficient or not, for computing averages of 
interest (the subject of this section); on the other, an efficient method for computing these averages (the subject of Taken 
together, these maneuvers produce a regularized vector field which one can evolve with giant strides using an appropriate 
numerical integrator. 

The calculations below largely foUowlffiil l|l882l ): they are repeated here partly to include the generalization to softening, 
partly to update Hill's derivation to modern notation and style, and partly to express the results in a vectorial, coordinate- free 
notation that is much more efficient for numerical integration algorithms. 

2.1 Regularized equations 

We consider a particle (called the 'planet', though it could be a star, moon, planetesimal, comet, etc.) bound to a massive 
central body (called the 'star', though it could be a black hole, planet, etc.). The planet follows a Keplerian orbit and is 
subject to relatively weak external forces, which cause time variations in the (osculating) Keplerian orbital elements. We 
express the equations governing these variations in terms of non-singular vectorial variables, then recover their secular average 
for the case of interest, namely when particles perturb each other's Keplerian orbits with their mutual gravitational attraction 
(Newtonian or softened). 

We denote by M* the mass of the star, and adopt the following notation for properties pertairung to the planet: 

• m mass of the planet, 

• r radius vector from the central mass, with r = \r\, 

• V = r velocity vector, 

• a the semi-major axis, 

• e the eccentricity, 

• n — ^G{M^ + m)/ the mean motion, with P = 2n /n the orbital period, 

• I, E, and / the mean, eccentric and true anomalies, respectively. 

When considering perturbations between planets, unprimed variables refer to the perturbing planet, and the primed counter- 
parts to the perturbed planet. 

In the planetary context, it is common to describe the perturbations to a Keplerian orbit usin g formulae that describe 
the variation of the orbital elements eccentricity, inclination, etc. (e.g., iMurrav fc DermottI Il999l ). These equations suffer 
from singularities at zero eccentricity and inclination, as well as at unit eccentricity. The singularities at zero eccentricity 
and inclination are traditionally handled by a canonical transformation to appropriate Cartesian coordinates for the orbital 
elements (the familiar trick for dealing with the singularity of cylindrical polar coordinates at the origin) . Such transformations 
work reasonably well in 'cold' settings where the eccentricity and inclination are small. Ne ar unit eccentri city (nearly radial 
orbits) the singularities can also be removed by an appropriate canonical transformation l|Tremainell200ll ). However, in hot 
environments such as the stellar cusp around a black hole the eccentricities and inclinations take on all values so one must 
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regularly switch between coordinate charts if one insists on working with elements such as eccentricity and inclination or 
even their Cartesian analogs. We initially implemented such switches, and found them costly and cumbersome. Instead, we 
chose to parametrize an orb it with the coor di nate- free (though redunda nt) angular momentum and eccentricity (or Laplace, 
or Runge-Lenz) vectors (see lHagihara|[l97ll or lMurrav fc Dermottlll999l ). 

We associate with an orbit three orthogonal unit vectors: z, along the rescaled angular momentum vector L = y/1 — z 
(normal to the orbital plane, scaled by mna^), x, pointing toward periapsis, and y — zXx, perpendicular to both. The 
eccentricity vector is 

. vxL , . 

A — r — ex. (1) 



The position vector of the planet is r = xx + yy where x = r cos / = a(cos E — e), and y = r sin / = aVl — sin E, with 
the radius r — a(l — ecosE). We introduce f, a radial unit vector, and the tangential vector t = zXr; explicitly 



(cos E — e)x + yl — e^sin E y 

r = cos f X + sin f y = 

1 — e cos E 

/ (2) 

— y 1 — sin E X + (cos E — e)y 

t — ~ sin f X + cos f y — . 

1 — e cos E 

We will also need the velocity vector v = xx + yy, where x = —na^ sin E/r and y = na^Vl — cos E/r (using E = na/r). 

The semi-major axis, angular momentum vector, and eccentricity vector of the disturbed (primed) planet respond to a 
perturbing acceleration /' on that planet as 

da' 2 



dt 

djn'a"' L') 
At 

dA' 

At 



r'xf 

The first equation expresses the change in the Keplerian energy resulting from the perturbing force, the second the change in 
angular momentum resulting from the associated torque, and with it the change in orientation of the orbital plane; the third 
governs the eccentricity and periapsis direction of the orbit, hence dictating shape and orientation in the orbital plane. 
We shall work with the radial (R), tangential (5*) and normal components (W) of the perturbing acceleration: 

R = r'-f, S = i'.f' W = z'-f. (4) 

Substitute into equation (|3]) to get 

da' _ 2(fle'sin E' + S y/ 1 - e''^) 
'At ^ n'(l - e' cosE') ' 

^^^^ = K,^^N^,y'^K,z', 

^ = A'^,x' + J^y,y' + K,z', (5) 

where 



N'^, = Wa'v/l-e'2 sin£', 
N'y, = -Wa'{cosE' -e), 



K, = ^^(l - e'cos£;'), (6) 
and 



^, _ yi - e'^[(4 cos E' ~ e' cos 2E' - Se')^ 2^1 - e'^ sin E' R] 
~ 2n'a'(l - e'cos E') ' 

[2(2 - e'2) sin E' ~ e! sin 2E'\ S - 2^1 - e'2 (cos E' - e') R 



y 



2n'a' (1 - e' cos E') 



K' = -^WsinE'. (7) 



n'a' 



These equations govern the secular evolution of an osculating Keplerian orbit in the vectorial representation. They are general 
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enough to accommodate a great variety of forces, conservative (external perturbers, such as sateUites, planets, stars, or 
black holes, on asteroids, planets, or stars) or not (drag, tides, dynamical friction), that may affect the motion of a particle 
whose motion is otherwise dominated by a massive central body. In this work, we are concerned with analytic averaging over 
softened mutual gravitational interactions in a near-Keplerian cluster of particles, to which we shall restrict our attention in 
the following. 



2.2 Averaged equations 

The softened gravitational force on the primed planet, m'/' — 

XX + yy + zz 



-V'"!", derives from the disturbing function 



<1> — —Gmm 



1 



(8) 



where Aj = •y/(x — x'Y + {y — y')'^ + (2 — z'Y + fe^, and h is the softening length (primed and unprimed refer to perturbed 
and perturbing respectively). The first term is Plummer's potential (which reduces to Newton's for & = 0) coupling the two 
particles, the second (the indirect term) is a non-inertial contribution resulting from pinning the reference frame to the central 
mass. Note that we soften only the interplanetary potential, not the potential between each planet and the central star. We are 
interested in the secular evolution of the osculating Keplerian orbital elements, which we obtain by averaging the equations 
of motion ((Tjl over the mean anomalies I and I' of the disturbing and disturbed planets. In other words we are after 



"da'' 






/•2tt 


/•2n 


da' 


At 


w 


47r2 


Jo J 


'0 


d^' 


AL'' 








p2TT 


AL' 


At 


w 


47J-2 


Jo J 


'0 


At 


AA'' 








n2TT 


AA' 


At 


w 


4^ 


Jo J 


'0 


At 



Al'Al. 



(9) 



Before proceeding further, we make a couple of simplifying remarks, (i) The indirect term in $ contributes a force r/r 
that averages to zero over Therefore we may drop the indirect term, (ii) Another classical result of secular dynamics is that 
Aa' /At averages to zero over /', since the average is proportional to the work of a conservative force over a closed path. Thus 
we may treat the semi-major axis a' and mean motion n' of the disturbed planet as constants of motion when examining 
conservative perturbations in secular dynamics (in the case of unsoftened gravitational forces this is Poisson's theorem, but 
the result requires only that the perturbing force is conservative and time-independent and hence holds for softened forces as 
well). 

Below, we shall demonstrate that the average of the direct force over I is analytic and expressible in terms of elliptic 
functions. For now, we denote by [R]i, [S]i, and [W]i the /-averaged radial, tangential and normal components of the direct 
force, and note that the i-averaged equations can be directly translated from their unaveraged counterparts in equations ((Gjl 
and (O, by replacing R, S and W by their respective I averages. 

Now consider the average over the orbit of the disturbed planet. For any function X{E'), 



X{E'){1 ~ e'cos E')AE'. 



(10) 



Inspecting equations ((Sjl and ^ we find that the average over I' will at most involve the second-order term of the Fourier 
expansion of [i?];, [S]i and \W]i in their periodic dependence on the eccentric anomaly E' , leading to: 



AL' 

At 

AA' 

At 



[KAw^' + [N'^,\^,v' + [K']urz'., 

[^^^,\,x' + [A',,\^,y' + [K']urz' 



(11) 



where 
and 



a'yJl-e''^{Wl ^\e'W^), 

-a'[{l + e'^)Wl-le'W!^-\e'W!] 

a'[{l + \e'^)S'',~2e'Sl + \e'^Sl\, 



(12) 



2n'a 



j^ml - e'Sl - 3e'S,") + 2^1 - e'm\\, 
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[A'y']w = ;^[2(2 - e'^)Sl - e'S'i - 2^1 - e'^Rl - e' R% 
" In a 

Here 

1 Z"^" 1 Z"^'" 

Rc=7^ [R]i COS nE' dE', R"s ^ TT \R]ismnE' dE' , (14) 

with similar definitions for S", ly", etc. 

For conservative forces, [da'/df];;/ = and this gives us the identity 



e R\ + \/\ — e'^Sc =0 if /' is conservative. (15) 

The cost of non-singular coverage of the full range in eccentricity and inclination using the vectors L' and A' has been to 
expand the phase space of an N-ring system from 4A'' (reduced from 6A'^ as a result of averaging over the mean anomaly and 
conservation of the semi-major axis a'), to the redundant 6A^-dimensional phase space spanned by the angular momentum and 
eccentricity vector of each ring; A'^ of the redundant components are removed by noting that L' and A' are orthogonal, while 
the other N reflect the relation between the magnitudes of the two vectors, L'^ + A'^ = 1. It is straightforward to show that 
[dL' -A' /dt]iii — and [d(L'^ + A'^) /dt]iii — are satisfied by equations (|li p - (|13|l if the perturbing forces are conservative. 



2.3 Analytic averaging 

We shall go over the crux of our exercise, namely the demonstration that Gauss's averaging works out perfectly for softened 
potentials. 

The average of the gravitational potential $ over the orbit of the perturbing (unprimed) planet can be written 

. , , , Gm f^" 1-ecosE 

) = -TT- / '^E , (16) 



27r 



where 



A? = - 2B cos{E - e) + C cos^ E, (17) 
with 

Ab — r'^ + + + 2aer' -x, 
Bcose = ar'-x + a^e, 
Bsine — ayA^-cP r' -y, 

C = ae^. (18) 
Note the inequalities 

Ab>a(l-e^)+h''>Q. (19) 

We pause briefiy to emphasize the seemingly innocuous absorption of the softening 6^ in Ai,. As far as Gauss's method 
is concerned, this is the only direct change that is brought about by softening. 
The average perturbing acceleration is 



■"6 



Gm f^" l~ecosE I . r ^ . ^ , 

/ dE -TT (r + aex — a\/ I — y smE — axcos E) 

27!" Jo ^6 



Gm f^^ , ^ 1 — e cos E 



„ , dE -= (Fo +Fisin£ + F2CosS), (20) 

27r /„ A? 



where 



Fo = —r' — aex, Fi = ay^l — fp- y, F2 = ax. (21) 
At this point, we write 

smE=—, cos£=— ; (22) 
introduce a new angular variable T (the perspective anomaly) defined by 
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sinT=^, cosT=^; (23) 

yo yo 

and demand that 

x = Qy (24) 

where Q is a 3 x 3 matrix with elements Qij, i,j — 0, 1, 2. 

The geometric motivation for this change of variable is described more fully in Appendix [X] Briefly, the orbit average [■]i 
can be thought of as a line integral around the orbit, with r and E the polar variables in the orbit plane for this line. We show 
in the Appendix that the orbit average is unchanged if we shift the line integral to a new contour, deflned by the intersection 
of any plane with the cone that contains the orbit and has its apex at the disturbed planet. The perspective anomaly is a 
polar coordinate in this new plane. The plane is later chosen to make the line integral as simple as possible. 

Since sin'^ E + cos^ E = sin^ T + cos^ T = 1 we must have 

xl-xl~ xl = 0, yl-yl-yl = 0. (25) 

The first of these simply states that x is confined to a cone; the second restricts Q to be a matrix that conserves the quadratic 
form x^ Mx, where x'^ is the transpose of the column vector x and M is the diagonal matrix diag (1, —1, —1). This requirement 
can be written 

Q'^MQ = M. (26) 

Thus Q is a pseudo-orthogonal matrix; in the language of special relativity, Q is a Lorentz transformation in two spatial 
dimensions, M is the Minkowski tensor, and the vectors x and y lie on the light cone passing through the origin. 
We define a diagonal matrix with complex coefficients, C = diag (1, i, i). Then 

C^ = (C*)'' = M, C*C = I, CMC = C*MC* = I, (27) 




where * denotes complex conjugation and I is the unit matrix. Let L = C*QC*. Then using equations (|26|l and p7l) we have 

L^L = C'Q^MQC* = C*MC* = I. (28) 

Thus L is orthogonal if Q is pseudo-orthogonal, and any pseudo-orthogonal Q can be written in the form 

Q = CLC (29) 

where L is orthogonal. Note that even though L may have complex coefficients, it is orthogonal {L~^^ — Lji) rather than 
Hermitian {L~^ = L*j). Because we allow complex coefficients the inner product of a vector with itself, x'^x — x-x — 5^^_q x1 
is not necessarily positive or even real. 

We may now re-write equation (|17p as 

xqA\ — AbXQ — 2BxoX\ sine — 2BxqX2 cose -|- Cx\ — x^Px, (30) 

where 

— Bsine — Bcose \ 

■ (-^1) 

C ) 

Gauss now asks for the pseudo-orthogonal matrix Q that diagonahzes the quadratic form for x'q/^1, that is, 

x'^Px = y^Q'^PQy = y'^Dy, (32) 

where D is some diagonal matrix. Using equation (|29p we re-write the quadratic form as 

y^CL^RLCy, (33) 
where 

CAb — iSsine — iScose \ 

-iSsine ° ' (3^) 

-iScose -C J 

Since R is symmetric, (i) its eigenvectors e-' can be chosen to be orthonormal, e-'-e'° — Sj^; (ii) the matrix L defined by 
Lij — e[ is orthogonal; (iii) L^'^RL = diag (Ai), where {Ai} are the eigenvalues of R. Thus equation (|33p can be written as 

y^CL^RLCy = y^Cdiag (A,) Cy = y^diag (Ao, -Ai, -A2)y, (35) 

which is the desired diagonalization (|32|l . with Q related to L by (|29p and D = diag (Ao, — Ai, — A2). 
The eigenvalues of R are the solutions of the cubic equation 



© 2008 RAS, MNRAS 000, 000-000 



8 J. Touma et al. 



y{\) = + (C- Ab)\' + (B^ - AiC)\ + B^Csin^ e = 0. 



(36) 



Using the relations p8|) it is straightforward to show that 

y{~C) = -B^Ccos^e<0, 
2/(0) = B^Csin^e>0, 

y{Ab) = ^^(Ab + Csin'e) > 0. (37) 

Furthermore from (|19p a^(l — e^) < ^i, so J/(A) must have three real roots, one between —C and 0, one between and 
(1 — e^) , and one between a^(l — e^) and A;,. It proves useful to label these roots so that Ao > Ai > A2; thus 



-C<A2<0<Ai< a\l - e^) < Ao < A. 
Relations between these roots include 

A0A1A2 = -B^Csin^e 
(Ao + C)(Ai + C)(A2 + C) 
A() + Al + A2 
AqAi + A1A2 + A2A0 
Explicit expressions are 

Ao = -2v/Qcos(i6»+ fyr) - i(C-^6) 



C 

-A.C. 



2VQcos(i6») - |(C- A) 



Al = -2VQcos(|e»- fyr) - i(C-^6 
A2 = 
where 

Q = i(C-Af,)'- i(B'- AC), 

i? = ^{C-Abf-\{C-Ai){B''-AiC) + ^B^Csm^e, 



27 
COS 



The eigenvectors of R are 



B sin e B cos e 



Ofc 1, 



At 



Afc+C 



where a^, which may be complex, is chosen so that e -e = 1 or 

B^ cos^ e 



-1 + ' 



+ ■ 



(Afc + C)2_ 

Using the relations (|39p this can be manipulated into a more compact form, 
Afc(Afc + C) 



(Afc — Ai)(Afc — Am) ' 
where I, m are the indices other than k from {0, 1, 2}. With the ordering (|38|l . 
Qq < 0, > 0, al > 0. 

We choose the phases of the eigenvectors so that 
ao = -i\/-ag, ai = --\/a?, 02 = 



Then 
Q = CLC 

/ 



^^- 



-OgB sm e 



/Ao 



QfjB sin e/Ai 



y^Bcose/(Ao + C) y^Bcose/(Ai +C) 



QjB sin e/A2 
a2Bcose/(A2 + C) 



(38) 



(39) 



(40) 



(41) 
(42) 

(43) 

(44) 

(45) 
(46) 



Ao(Ao+C) 



(Ao-Ai){Ao-A2) 

rB sin e 



Ai(Ai+C) 



Aq+C 



Ao(Ao— Ai)(Ao — A2) 
Aq 

(Ao+C)(Ao-Ai)(Ao-A2) 



(Ao-Ai){Ai-A2) 

B sin e 



|A2l(A2 + C) 
(Aq-A2)(Ai-A2) 



Ai+C 



Ai(Ao-Ai){Ai-A2)-' 



A2+C 



A2l(Ao-A2){Ai-A2) 



B cos e 



Al 



(Ai+C)(Ao-Ai)(Ai-A2) 



B cos e 



IA2 



B sin e 



rB cos e 



(47) 



(A2 + C){Ao-A2){Ai-A2)' 
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The determinant of this matrix is 

det (Q) = det (L) = det (e^) = s where s = sgn [sin(2e)] = ±1. (48) 

We now examine the relation between the eccentric anomaly E and the perspective anomaly T in more detail. From the 
(0, 0) component of the equations Q"'"MQ = M and QMQ""" = M we have 

Qlo + Qlo = Qoi + Qo2 = Qoo - 1- (49) 
Set 

^ = V Qoo - 1' sm(t>=-^, cos(f>=—, smi/)=— , cosi/' = (50) 
it is straightforward to show that Qoo > 1 so A is real and positive. Then from equations (I22p and (|23p 
Axo + Qioxi + Q20X2 = Axo[l + cos{E - (/>)], Ayo + Qoiyi + Qo2y2 = Ayo[l + cos(T - V)]- (51) 
From equation (I24|l we also have 

Axo + Qioxi + Q20X2 = yo(AQoo + Qio + Q20) + yi{AQoi + QioQii + Q20Q21) + y2{AQo2 + Q10Q12 + Q2oQ22)- (52) 
Using the (0, 0), (0, 1) and (0, 2) components of the equations Q^'^MQ = M we can simplify this to 

Axo + Qioa;i + Q20X2 = (Qoo + A)(Aj/o + Qoiyi + Qo2y2), (53) 
so from equations (|51|l we have 

a;o[l + cos(£ - 4>)\ = yoiQoo + A)[l + cos(r - V)]- (54) 
Similarly, 

Q20X1 - Q10X2 = Axo sm{E - (t>), Qo2yi - Q01J/2 = Ayo sin(r - 1/)), (55) 
and 

Q20X1 — Q10X2 = yi{Q2oQii — Q10Q21) + y2{Q2oQi2 — Q10Q22). (56) 
Using equations (|39|l and (|47p it is straightforward to show that 

Q20Q11 — Q10Q21 ~ sQo2, Q20Q12 — Q10Q22 ~ —sQoi (57) 
where s = ±1 is defined by equation (|48|) . Thus 

Q20X1 — Q10X2 = s{Qo2yi — Qoiy2), (58) 
so 

xosm{E — (f>) = syosin{T — Tp). (59) 
Combining this result with (I54|l . 

tani(r-V) = s(Qoo + A)tani(£'-0) = s (^Qoo + ^Qoo - l) tan i(£ - </>). (60) 

This result relates the perspective and eccentric anomalies and shows that the perspective anomaly circulates through 27r (if 
s = 1) or ~2n (if s = —1) when the eccentric anomaly circulates through 2n. Moreover the differentials are related by 

dT = .(Qoo + A) --;f;^-;i di. = ^l^^di. = s^Ji.E, (61) 
cos2 i(_E - 0) sm[E - 4>) yo 

where the last equality employs (|59|l . 

We may now combine equations (|20p . (|22p . (|23p . (|30p . and (|35p into an expression for the orbit-averaged perturbing 
acceleration at position r: 

f2n 



lf']i = ^ f^E ^ ^^3°^ ^ (Fo + f 1 sin g + f 2 cos g) 

Gm /"^"^^ 2;o(xo - ex2)(Foa;o + f 1X1 + ^2x2) 



27r (ylijXQ — 25x0X1 sin e — 25x0X2 cose + Cx|)^/2 

Gm Z"^" , yo Y.%oiQo] - eQ2j)yj k=o 2 ^jQjkyk 



27r X [Ao-Ai(yi/j/o)2-A2(y2/j/o)2]3/2 ' 
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Replacing yi/yo by sinT and 2/2/2/0 by cosT, and dropping terms that integrate to zero, we have 

(63) 



(64) 



[f']l = 


Gm /"""^y Fu + Fvsin^T 
27T X [Ao - A2 - (Al - A2) sin^ r]3/2 


wh6r6 








2 


2 


Fu = 


E- 










the F 


j are 


defined by equation (|2H). and 


Uo = 


: Q 


'00 ^ e(5ooQ20 + Q02 ~ eQo2Q22 




= Q 


'ooQio — &Q10Q20 + Q02Q12 — eQi2Q22 


U2 = 


: Q 


'00(320 — eQio + Q02Q22 — e(522 


Vo = 


: Q 


*oi ~ eQoiQ2i — Q02 + eQo2Q22 


Vi = 


-= Q 


'oiQii — eQiiQ2i ~ Q02Q12 + eQi2Q22 


V2 = 


-= Q 


'01(521 — ('Q21 ^ Q02Q22 + e(322- 



(65) 

The integrals can be expressed in terms of the complete elliptic integrals E{k) = J^^^ di9(l — fc'^ sin'^ "^Y^^, K(k) — J^^'^ di9(l — 
P sin2,9)-i/2. 

'"'^^ dT _ E{k) /""/^ sin^TdT _ E{k) K{k) ^^^^ 



[l^k^sm^Tf^ l-fe^' 7o (l-fc2sin2r)'/^ k^{l-k^) P ' 
Thus 

'^']' = ^7r^f °^ [(fcV. + Fv)i?(fc)-(l-fc^)F.j^(fc)], fc' = ^^. (67) 

TT (Aq — Aij(Ai — A2j Ao — A2 

The averaged radial, tangential and normal components of the force ([-R];, [S]i, \W]i) are found by replacing /' in equations 
Q by its averaged value (|67p . 

Thus, one of the averages is done, the force of a softened ring at a point is calculated. To get the fully averaged secular 
equations (|12|l . a second average over the perturbed ring itself is required, as given by equation p4|) . The second average is 
accomplished by numerical quadrature as described in the following section. 



3 REMARKS CONCERNING IMPLEMENTATION 

To be sure, we are not yet at the point where we have an optimal N-ring code, partly because the improvements afforded by 
averaging lessen the urgency to optimize, partly because we prefer to explore novel dynamical phenomena rather than novel 
algorithms. Obviously, a careful optimization will eventually be necessary, but for now we simply outline the steps involved in 
the efficient implementation of the Gaussian ring algorithm, emphasizing difficulties encountered, and possible ways around 
them. We are engaged in optimizing and parallelizing the code, with a view to simulations of stellar clusters surrounding 
massive black holes. 



3.1 Numerical averaging 

In order to perform the numerical average over the orbit of the perturbed planet, which gives us equations pi|l - p3| l for the 
evolution of its angular momentum and eccentricity vectors, we need to compute the first three Fourier coefficients in the 
eccentric anomaly E' of the force vector (i?!?, Rl, . . . , ; see eg. I14|l . Direct summation at points equally spaced in E' (i.e., 
a discrete Fourier transform) will do in most cases, and we have adopted this approach in the calculations presented below. 

The required number of points K for the evaluation of these Fourier coefficients depends both on the error one is willing 
to tolerate and on the eccentricities and separations of the rings. The Euler-Maclaurin summation formula tells us that the 
accuracy of the evaluation will improve with increasing K faster than any power of 1/K, typically as exp{—aK) for some 
constant a > 0. Of course, during close approaches between rings this desirable asymptotic behavior may only occur at rather 
larg e K, and th e constant a may be disappointingly small. 

Imi] (|i882l ). who was mainly concerned with planetary motion (nearly circular, nearly coplanar, well-separated rings), 
was happy with K — 12, even 8. More generally, he pointed out that 'if the number of these values be even, the order of the 
error. . .will be the same as that of the power of the eccentricities or mutual inclination of orbits, whose exponent is one less 
than the number of these values'. We can re-state Hill's conclusion-which we have confirmed using expansions of the disturbing 
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functions in eccentricity and inclination-more precisely as follows: Consider two rings with eccentricities and semi-major axes 
Cj, Uj, j = 1, 2, and mutual inclination /. Let 5 — \a2 — ai\/a where a — ^(ai + 02), and let e — max(ei, 62). Then if e <C 5, 
7 <C 5, the fractional error in the time derivatives of the eccentricity and angular-momentum vectors is 0[max{e/S, I /S)^~^] 
as e.I —> This result also implies that the implementation of Gauss's method using numerical quadrature with K equally 
spaced points is formally equivalent to a secular perturbation theory in which the Hamiltonian is expa nded to first order in 
the planetary masses and K^^ order in the eccentricities and inclinations (e.g., iLiiaert fc Henrardir2008l . for K — 12). 

Going by this estimate, Hill's use of K = 12 would yield a relative error of 10"^" in the time derivatives of the eccentricity 
and angular- momentum vectors of a two-planet system consisting of Jupiter and Mercury (eccentricity 0.2); numerical experi- 
ments on this system, with varying orientations, yielded relative errors in the range 10~^-10~'^^ in reasonable agreement with 
this estimate. To achieve a similar accuracy for the system consisting of Jupiter and Halley's comet our experiments indicate 
that one needs K ~ 200-400, and this again depending on argument of node and perihelion, while keeping the other elements 
fixed. A softening length b = O.Ola or fe = 0.1a can cut this number by about a factor of two or four, respectively. The stellar 
systems which we will be dealing with have higher eccentricities and inclinations than anything Gauss or Hill envisioned, as 
well as frequent ring crossings. One can get away with K = 16-32 for nearly circular rings, while achieving a relative error of 
10"^" or smaller in the evolution equations. For nearly radial orbits, we found that we needed up to if = 512 for the same 
level of accuracy. 

It is useful to have a criterion for the accuracy of the numerical average, which can be used to provide an adaptive choice 
for K for each ring pair. Hill's criterion works when e 6, I <^ 5, but not for near-radial or crossing rings. We devised an 
alternative criterion which is equivalent to Hill's in the limit where Hill's works. We rely on the observation that, had the 
averages been exact, the equation governing the evolution of a ring's semi-major axis would be [da'/dt];j' = (see discussion 
following eq.O, which in turn implies the identity p5|l . Therefore we may simply choose K adaptively so that the numerical 
evaluation of the absolute value of the left side of (I15|l . divided by n'^a to make a dimensionless quantity, is less than some 
specified tolerance tquad ('quad' for quadrature) for each ring pair. We start with K — 16 and proceed through doublings. 
Note that Richardsonian extrapolation is not useful here, because the asymptotic convergence is faster than any power of 
1/K. What is an appropriate tolerance? A tolerance tquad translates into a tolerance on the relative error in semi-major axes 
per secular or precession period that is no worse than (A*' — l)equad. This in turn translates into into a relative error over M 
secular periods no worse than AMequad. Presuming one is interested in simulating a system of 10"" particles over 100 precession 
periods, and one can live with cumulative relative error in the dominant Keplerian energy of lO"*, then 10"^^ should 

be amply sufficient for such a bound. For most ring-ring interactions, such a tolerance is achieved with 64 to 128 sectors per 
ring at a softening of 0.01 times the typical size of a ring. There are of course likely to be some ring pairs for which sampling 
of 16 to 32 times these values may be necessary, and in these cases adaptive sampling proves essential. 

So far we have described how to carry out numerical averaging over the perturbed ring using equally spaced points in 
eccentric anomaly. However, we are aware that more sophisticated quadrature routines may be better, particularly at ring 
crossings when the softening length b is small. We have conducted limited experiments with the QAGS adaptive integrator 
from GSL (intended to handle integrable singularities by adaptive application of the Gauss-Kronrod [Gauss again] rule). In 
the rather extreme case of two coplanar, intersecting rings, with a softening of only 0.001 times the orbital radius, QAGS 
required 1500 judiciously placed integration points to achieve a relative error of 10~^^; ten times more points were required by 
uniform sampling to achieve the same tolerance. The factor of 10 drops to 2 for a softening length which is ten times larger. 

3.2 Numerical integration 

The rings are evolved in time using a Bulirsch-Stoer (BS) integrator. The BS integrator is not tailored to this problem, but 
for a reliable, brute force exploration, it is efficient and simple to implement. The ideal of course would be to exploit the 
Lie-Poisson structure of the equations to construct a geometric integrator with an adaptive timestep that would handle the 
broad range of dynamical time-scales in typical stellar and planetary systems. This is certainly a noble undertaking, which 
we put off to a later stage, and we will be more than happy if others beat us to it. In any case, in our BS integrations, we set 
a tolerance tint for the error in the integrated vector which is similar to the tolerance tquad that we demand in the numerical 
quadratures. In most cases, we are quite happy with a tolerance of eint ~ 10~* (at which we live with step sizes on the order 
of 0.1 secular periods); we do push things further in the simulations below, and evolve some configurations with a tolerance 
of tint = lO"'^^, and consequently smaller timesteps on the order of 0.01 secular periods. 

Highly eccentric orbits are likely to be common in ma ny nearly Keplerian syste ms (e.g., stellar clusters around black holes, 
cometary clouds). In the plane one has the lens orbits of Sridhar & Touma (1999); related families exist in three dimensions 
and were studied in oblate and prolate potentials bv lSambhus fc Sridhar (2000 ). The averaging process we describe here is 
well-behaved for nearly radial orbits (e — > 1), as are the equations of motion for the eccentricity and angular-momentum 
vectors-which is one reason why we chose to work with vectors rather than orbital elements, which usually give singular 

* Hill limits this result to the case where K is even. In fact, for rings with small but non-zero eccentricity the result also holds for odd 
K. Hill's qualification apparently relates to the case of inclined, circular rings: in this case odd K yields much faster convergence than 
even K, particularly if one point lies on the mutual node. 
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equations of motion as e — > 1 (|Tremaindl200ll ). Nevertheless we have found that our BS integrator struggles with nearly radial 
configurations, slowing the integration by a factor of 10 to 1000 in some extreme cases. 

Part of the slow down is dynamical, and will be faced by any adaptive integration scheme. The other, and more critical, 
part is a consequence of the violation of the constraint + = 1 in the course of the numerical evolution; it manifests as a 
square-root singularity in the equations of motion (|ll|| - (|13p . and forces a slowing down in the integration, as steps are reduced 
to avoid that singularity. One can imagine a number of solutions to this problem. As suggested above, the best would be to 
design a geometric integrator which is adapted to the Lie-Poisson structure of the equations at hand, and hence preserves the 
integral in question automatically. Less elegant, but still effective, would be to enforce the constraint + — 1 hy hand, 
renormalizing one or the other vector at regular intervals to eliminate the small errors introduced by the integration scheme. 
Without modifying the BS integrator, one could impose a reflecting boundary condition at high eccentricity. There are also 
physical processes that mitigate the problem in many systems: these include general-relativistic precession (Appendix|B)) and 
tidal disruption, collisions, or swallowing of stars by the central object. A careful examination of the proper implementation 
of these and related alternatives would take us too far afield. In the ring simulations described below, which include general- 
relativistic precession, we encountered an occasional near-radial ring that slowed down the integration of the whole cluster. We 
found that this problem could be removed either by brute-force removal of highly eccentric orbits (e > 0.999) or by exploiting 
the redundancy in the present vectorial formulation. In the second solution, we simply used the eccentricity vector A to 
solve for the eccentricity at high eccentricities and the angular-momentum vector L to solve for the eccentricity otherwise. 
The resulting algorithm still slowed down at high eccentricity passages, but only by that dynamical factor which reflects the 
higher precession rates at these eccentricities, boiling down to a factor of 10 decrease in performance during these phases. We 
found this intermittent slowing down-about a factor of ten when a ring (or a population of rings) developed mean eccentricity 
greater than 0.8-tolerable in the cases we studied, and believe it manageable with larger systems, when the code is eventually 
parallelized. 

3.3 Evaluation of elliptic integrals 

To compute the elliptic integrals one could either iterate using the arithmetic-geometric mean (Gauss once again), interpolate 
in a table, or rely on Chebyshev polynomial expansions. We tried all three methods and found that smart searches through 
tables were twice as costly as Chebyshev expansions, which were marginally better than the arithmetic-geometric mean (for 
an absolute error goal of 10"^''). In our current implementation, elliptic integrals are calculated with the arithmetic-geometric 
mean. 

3.4 Additional physics 

Nearly Keplerian stellar systems are often subject to additional effects with substantial secular contributions. We can easily 
accommodate the average contributions of tides (Galactic for cometary clouds, bulge/halo for stellar clusters around a central 
black hole). We can also account for dynamical friction contribut ions, in a manne r analogous to related efi'ects which one of us 
has handled in studies of the history of the lunar orbit ifTouma fc Wisdomlll998l ). Finally, we can also account for the secular 
effects of general-relativistic corrections f Appendix [B|) . 



4 THE ALGORITHM TESTED 

In this section, we test the Gaussian ring algorithm by comparing its results to the behavior of the equivalent unaveraged 
N-body system, integrated with the same Bulirsch-Stoer (BS) algorithm at equivalent tolerance. Alternatively, one may match 
numerically averaged solutions with results of analytical averaging, an exercise which we leave to a forthcoming discussion of 
analytical theories of softened systems of coupled rings (Mroue & Touma, in preparation). 

The desired accuracy of the calculation is specified by two parameters: the first is equad, the fractional accuracy of the 
numerical averaging, which was discussed in ^ and the second is eint, the tolerance specified in the numerical integrator for 
the ordinary differential equations. 

A technical point is that when comparing secular computations to unaveraged calculations, one should choose the 
initial conditions of th e unaveraged system to ensure that the comparison is not affected by high-frequency oscillations 
(ILaskar fc Simonlll988l ). This can be done by experimenting with various initial conditions for the unaveraged calculations 
to find one at which the high-frequ ency oscillation is near a node, or by techniques such as 'warming up' the unaveraged 
integration (|Saha fc Tremainelll992l ). In the present paper, however, we simply use the same initial conditions for both N-body 
and averaged calculations. 

The simulations below involve stars orbiting a 10^ M© black hole. Distances are measured in parsecs, time in years, 
and mass in units of the black-hole mass. In these units, Kepler's law relating the period P and semi-major axis a is 
P = 2.96 X 10* yr (a/pc)^^^. Orbital inclinations are measured with respect to the x-y plane of a right-handed reference 
frame, with origin at the central black hole; longitudes of node and periapsis are referred to the a;-axis of this frame. The 
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dominant relativistic correction (apsidal precession due to the central mass) is sometimes included, as described in Appendix 

m 



4.1 Softened Kozai oscillations 

To start out, consider a system that represents a star orbiting one of the components of a binary black hole at the centre 
of a galaxy. The black-hole companion is represented by a fixed ring having mass Mo — 1, semi-major axis Oo — 10, and 
eccentricity eo = 0.5 (recall that semi-major axis is in units of parsecs and mass in units of the mass of the primary black 
hole, lO'^ Mq). The fixed orbital plane of the black hole binary serves as the x-y plane of the reference frame, with the a;-axis 
taken perpendicular to the line of apsides, the j/-axis pointing to the periapsis, and the z-axis along the orbit normal. The 
star is represented by a ring with Mi = lO"'^, at = 0.1, and initial eccentricity — 0.01. Its orbital plane is initially inclined 
by li = 60°, cuts the plane of the binary along the a;-axis, and has a periapsis which is 90° from the node. The softening 
is O.lfli; relativistic precession is switched off. We first carried out a direct (restricted) three-body integration (BS tolerance 
tint ~ 10~^^) with these initial conditions, supplemented by randomly chosen mean anomalies. Then we performed a number 
of ring simulations for comparison. At a fixed sampling of 16 sectors per ring, a maximum error of 10~^^ in the Keplerian 
energy was comfortably achieved. The Gaussian ring algorithm follows the secular evolution with giant strides ranging from 
1.4 X lO'' (tint = 10"^*) to 4.5 X lO'^ yr (tint = W'^)- 

Below, we compare results of the secular ring computation at eint ~ 10~^^ with the benchmark particle simulation. In 
this rather conservative run, the timestep in the averaged simulation (~ 1.8 x 10^ yr, 1/40 of the Kozai period) was 2 x 10''' 
times larger than in the direct three-body integration (~ 90 yr, 1/10 of the orbital period). 

Neither energy nor angular momentum is preserved in the restricted, eccentric, three-body calculations. In the averaged 
dynamics, the star's energy is conserved but its angular momentum is not. The energy conservation provides a check on the 
accuracy of the Gaussian ring algorithm: the ring's secular energy was conserved to a maximum fractional error of 4 x 10"^", 
over a hundred Kozai cycleqf|. 

The orbital elements of the star from the three-body and secular codes are shown superimposed in Fig. [T] The eccentricity 
and inclination variations are due to Kozai oscillations. To the naked eye, the runs are practically indistinguishable in this 
figure. A closer look reveals two types of (minor) differences between the two codes, as shown in Fig. (2] (i) small-amplitude 
variations over the binary's orbital period, which average out in the secular dynamics; (ii) relatively larger variations with 
the same period as the Kozai cycle. The latter are mainly due to differences in the initial conditions between the three-body 
and secular codes; we were able to reduce them significantly by using averaged initial conditions for the secular dynamics as 
described above. While the differences in initial conditions are caused by second-order terms in the masses, the time evolution 
of differences is not; in fact, ring simulations with slight variations in initial conditions revealed similar (qualitative and 
quantitative) time variations in the differences of their orbital elements. 



4.2 Counter-rotating rings 

The lin e arized dynamics of softened, planar, counter-rotating, nearly Keplerian systems of stars and gas was considered by 
iToumal (|2002l ). There it was shown that (i) such configurations are prone to violent m = 1 (lopsided) instabilities; (ii) the 
bifurcation to instability, the growth rate and the pattern speed are functions of softening; (iii) the instability is associated 
with interaction of modes of opposit e energy. 

The secular theory employed bv lTounia l|2002l ) was based on a Hamiltonian that was truncated at O(e^). Mroue & Touma 



(in preparation) extend the analysis to 0(e ), and provide a description of the unfolding of the bifurcation, the appearance of 
uniformly precessing ring configurations, and the secular dynamics around them. The Gaussian ring algorithm offers an ideal 
tool for exploring this secular dynamics to all orders in inclination and eccentricity. Results on the planar two-ring problem 
will be presented in Mroue & Touma. Here, we further test our method by following the development of a nearly coplanar, 
counter-rotating system of two rings. 

Specifically, we follow two equal-mass rings, in initially nearly circular (0.01 eccentricity), nearly coplanar (1° inclination), 
counter-rotating orbits. We consider a case in which Mo — Mi — 10~®, flo = 1 and Ui — 0.5, with softening b — 0.3ao and 
relativistic precession switched off. This is a linearly unstable configuration, in which the eccentricity of both rings is expected 
to grow (note that reducing the softening length below 0.26ao switches off the instability, at this ratio of semi-major axes). 
The evolution is shown in Fig. [S] where one observes oscillations in the eccentricity from near-circular to near-radial and back. 
Even more dramatic is the inclination dynamics: the rings experience repeated flips between prograde and retrograde near 
the peaks of the eccentricity oscillations. The cause of this out-of-plane instability is not fully diagnosed, though it smells 
like Kozai cycles. To be sure, this behavior, which may stimulate bending instabilities in near-Kepleriari discs, is markedly 



° By secular energy, we mean the value of the averaged Hamiltonian that generates the secular dynamics. Relative errors in the total 
energy come out better (by several orders of magnitude) than one might expect for the given value of ei^t and tqu^d ; the reason is that 
the denominator is then dominated by the Keplerian energy —^GMi, /a, which is irrelevant to the actual integration. 
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Figure 1. Osculating elements of a star in a black- hole binary, as determined from a direct three-body integration (solid curve) and 
from the Gaussian ring algorithm (crosses) over a single Kozai cycle. 



different from the behavior of counter-rotating configurations studied in the context of the three-body problem. There, it 
has been understood for decades that counter-rotating configurations are rather stable, even more so than their prograde 
counterparts. It appears that softening (or heat) modifies the secular frequency spectrum to a point that angular momentum 
exchange between rings leads to an instability, which first promotes eccentricity growth to a near-radial configuration, then 
inclination growth to fiipped orientations, and back. 

In the course of these acrobatics, the secular algorithm satisfied an averaging tolerance tquad ~ 10~^^, by working its 
way up to 128 subdivisions per ring (from a minimum of 16) in the high eccentricity /inclination phases of the evolution. The 
dynamics was followed with a mean timestep of about 1/100 of the oscillation period, dictated by an integrator tolerance 



Cint 



10" 



At this tolerance, the relative error in the secular energy was no more than a few times 10 over 100 oscillations; 



the fractional error in the total angular momentum was bounded by 3 x 10~^^ over that same period. 

The reader may wonder how well we do against particle integrations in these vigorously unstable configurations. The 
situation is illustrated in Fig. |4] there is good agreement over a cycle, giving way to quantitative disagreement but similar 
qualitative trends. This discrepancy has little to do with proper adjustment of initial conditions. Rather, second-order effects 
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Figure 2. Difference between the osculating elements of a star in a black- hole binary, as determined from a direct three- body integration 
and from the Gaussian ring algorithm, over ten Kozai cycles. 

in the masses, which are neglected in ring simulations, are amplified in close encounters that result in the large eccentricity and 
inclination phases of this system. This is particularly apparent in Fig.O where jumps in particle semi-major axis accumulate 
to a few percent in less than five cycles. These jumps, which occur at close encounters/crossings between particles, are not 
accounted for in the secular representation of the dynamics. They bring about shifts in the secular frequencies which are 
responsible for the discrepancies observed in Fig. ijfl . 



In comparisons between ring and particle simulations it should be borne in mind that there are several possible particle analogs to 
a ring. Here, we associate a single particle with a ring, to highlight the discrepancies that the practitioner should keep in mind when 
studying secular few-body dynamics, say in the context of planetary systems. Alternatively, one could think of a ring as a guiding-center 
trajectory that provides the mean phase-space position of a distribution of particles, with total mass equal to the mass of the ring and 
spatial extent set by the ring's softening; this point of view is more natural to stellar clusters in galactic nuclei. Finally, the ring could 
represent a point in the five-dimensional space of actions and resonant angles, i.e., a component of a phase-space distribution function 
that is uniformly distributed in mean anomaly. The error associated with Gauss's method depends on the analogy being used. 
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Figure 3. Evolution of a counter-rotating system of two nearly circular, nearly coplanar rings. Top: eccentricity oscillations of the outer, 
initially prograde, ring (left), and of the inner, initially retrograde, ring (right). Bottom left: variations in the inclination of the outer 
ring. Bottom right: an expanded view of the exchange of angular momentum and inclination between the inner (solid) and outer (dashed) 
ring. 



Finally, while softening/heat can induce instability in a counter-rotating system of rings, general-relativistic precession 
can suppress the instability. The importance of normally negligible relativistic effects is already appreciated in studies of Kozai 
oscillations, (e.g., Holman et al. 1997), and we demonstrate that the same applies in this case by switching on relativistic 
precession in the same configuration (Fig. [Gjl . 



5 COUNTER-ROTATING SYSTEMS OF RINGS 

We conclude our battery of tests with a modest, yet challenging, look at systems of counter-rotating rings. For some time now, 
it has been re c ogniz ed that counter-rotation is likely to be associated with lopsidedness in self-gravitating systems of gas and 
stars. iToumal l|2002l ) pointed out that counter-rotating instabilities are at work in nearly Keplerian systems and speculated 
that non-linear evolution of unstable discs may lead to the formation of eccentric discs such as the one observed in the 'double' 
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Figure 4. Eccentricity and inclination of inner ring from three-body integration (dashed) and averaged integration (solid). Relativistic 
precession is not included. 
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Figure 5. Jumps in the semi-major axis of particles in the three-body integration during close encounters, in the course of oscillations 
caused by the counter-rotating instability. 

nucleus of M31 (|Tremainelll995l : iPeiris fc Tremainell2003l ). We limit our discussion to a couple of runs, which are illustrative 
of both the challenges faced by the algorithm, and the typical dynamical fate of such systems. 

We follow the dynamics of a system of 100 rings, half of them prograde with a total mass of 10"'' relative to the central 
point mass, and the other half retrograde, with the same semi-major axis distribution as the prograde rings, and 20 per cent 
of their mass. In this illustrative example, we picked the rings from uniform distributions in semi-major axis (over the interval 
[1, 1.2] pc), eccentricity (over [0, 0.3]), inclination ([0, 10°] for prograde and [170°, 180°] for retrograde), argument of node and 
periapsis (over [0, 27r]). The softening length was b — 0.1 pc. 

The integrator used a tolerance tint = 10~*, which gave a mean timestep of 10* yr; the timestep reached a maximum of 
6 X 10* yr and dipped as low as 100 yr for brief periods. The numerical averaging was performed with a minimum of if = 128 
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Figure 6. The stabilizing effect of general-relativistic precession as evident in the tame behavior of eccentricity and incUnation of the 
same ring system. 



points per ring and a maximum ol K = 512, enforcing an averaging error tquad ~ 10"^". The system was followed for 30 Myr. 
During this period, the secular energy was conserved to within a fractional error of 10~®, with ten times smaller error in the 
total angular momentum. The run took about two weeks on a single 3.2 Ghz processor. An order of magnitude improvement 
in performance can be achieved if one is wiUing to live with a tolerance eint = 10~®. 

The evolution of the system in eccentricity and inclination is shown in Fig. [7] In the first 5 Myr the system undergoes 
significant growth in eccentricity, at near constant inclination; by 5 Myr the mean eccentricity of the prograde population 
is around 0.6 and the mean eccentricity of the less massive retrogra de population is larger at 0.85. This result is not new: 
a closely related planar system is known to be linearly unstable (|Toumal 120021 '). with the lighter/retrograde component 
experiencing a more significant growth of eccentricity than the heavier/prograde component. It is the later fully non-linear 
behavior that we can confidently follow with the Gaussian ring algorithm for the first time. Following the initial growth of 
eccentricity, and in a manner reminiscent of the behavior witnessed in the previous section for two rings, the rings in the 
retrograde/lighter component experience a dramatic growth and spread in inclination, practically reversing orientation by 
7.8 Myr. The mean eccentricity of both populations increases slightly during this fiipover phase, with the retrograde rings 
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Figure 7. The instability viewed in the (e,i) plane. Prograde rings are represented by filled circles, while retrograde rings are open 
circles. As the time increases, from top left to bottom right, one witnesses the growth of ecceiitricit j' at near constant inclination, followed 
by the flip-over of the retrograde rings, then their attempted return to their initially coplanar configuration. The more massive prograde 
rings experience a significant net growth in their mean eccentricity, while remaining clustered close to their initial inclination. 
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Figure 8. Evolution of tlie projected and smoothed surface density of the full system. The projection is on the reference x-y plane, 
lumping together contributions from prograde and retrograde orbits. The central point mass is at the origin, marked by a cross. Color 
refers to density levels. An m = 1 mode grows as a result of the counter-rotation instability, capturing rings in its way, with the prograde 
population maintaining the precessing single-lobed pattern, all the while the initially retrograde population prepares for disintegration. 
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Figure 9. The projected surface density evolution of tlie rings over longer times, when viewed in negative color. The highly eccentric and 
dispersed population of initially retrograde rings projects onto a lightly hued halo, around the heavier lopsided configuration of prograde 
rings, captured in a precessing m = 1 mode. 



reaching a mean of 0.9 around 6 Myr. Then, the initially retrograde rings try to regain their home-base; they almost succeed 
around 10.3 Myr, but suffer significant losses along the way. At this point their mean inclination is about 140°, and their 
mean eccentricity has decreased to 0.6; the prograde rings are still in low-inclination orbits, with a mean eccentricity of 0.45. 
The prograde rings are also bunched in the conjugate angles, the longitude of periapsis in particular. This is apparent in 
the associated spatial densities, shown in Fig. [S] The non-linear saturation of the instability leads to an oscillating, eccentric 
ring composed of prograde particles, with a mean inclination of 7.5°, a mean eccentricity of around 0.6, and an increasingly 
steady prograde precession with a frequency 0.8 rad Myr~^ (Fig- IS]). Meanwhile the culprits behind the instability, namely the 
initially retrograde rings, are dispersed into a population with a wide spread in inclination, node and argument of periapsis, 
and a mean eccentricity of 0.7 (Fig. [TDJ. 

In short, the mode, which is ignited by the instability, grows through resonant capture of both prograde and retrograde 
populations, then saturates as the retrograde population (which has a lower angular momentum) escapes from resonance, and 
disperses through (partially stochastic) phase mixing; this leaves a captured prograde population which maintains a lopsided 
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Figure 10. The dispersal of the retrograde rings (open circles), and the increased cohesion of the prograde rings (filled circles), as seen 
in the (e, i) plane. 



m = 1 mode that is increasingly steady in its (prograde) precession. Just for fun, and once the mode was fully established, 
we removed all initially retrograde rings, and left the prograde population to itself. Perhaps not surprisingly, the captured 
prograde population, shown in Fig. 1111 maintained its cohesive state, now precessing uniformly with a slightly higher frequency 
of 0.9 rad Myr~^ for at least four precession periods. At a mean timestep of 1.1 x 10* yr, the integration kept fractional errors 
in energy and angular moment um below 10~^. On e coul d think of the resulting mode as an analog of the lopsided mode 
which was artificially excited bv I Jacobs fc SellwoodI (|200ll ). The main differences are: (i) the Jacobs-Sellwood simulation was 
a two-dimensional N-body simulation, while ours is a three-dimensional N-ring simulation; (ii) the Jacobs-Sellwood simulation 
used A'^ =100,000 while ours used A'^ — 100 (each of our rings should be regarded as representing many stars); (iii) the Jacobs- 
Sellwood simulation formed an eccentric disc using eccentric initial conditions, while we relied on counter-rotation to generate 
an eccentric equilibrium through the non-linear evolution of an instability. 

A careful treatment of the growth and saturation of the counter-rotating instability would take us too far afield. It is in 
fact the subject of an upcoming study by Touma and Kazandjian (in preparation), who perform sufficiently resolved N-body 
simulations to permit a detailed phase-space dissection of the dynamics in question. Their results, which are based on long 
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Figure 11. The evolution of the prograde rings after the retrograde rings were removed, leaving a slightly thickened, lopsided, and 
uniformly precessing annulus. 



term N-body simulations of up to 500,000 counter-rotating particles around a dominant central body, are in remarkable 
qualitative agreement with the salient features of the 100 ring simulations reported here. Still, the full fledged confrontation 
of our algorithm with a conventional N-body code will have to wait for ring simulations that are sufficiently resolved (10^ 
rings, pursued over a hundred precession periods) to permit detailed phase-space analysis, and comparison of the dynamics. 
These simulations are surely within the reach of a parallel version of our algorithm, running on a state-of-the-art cluster, and 
we hope to report on such undertakings in the near future. 



6 DISCUSSION 

The algorithm we have described computes the secular or orbit-averaged evolution of systems of particles in a Keplerian 
force field and interacting via softened gravity. It is accurate to first order in the masses of the particles, and to all orders in 
eccentricity and inclination. The order for a system of A'^ particles, which is 6A'' in a conventional A^-body code (assuming that 
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momentum conservation is used to solve for the motion of the central body), is reduced by 27V as a result of averaging and the 
consequent conservation of semi-major axis. However, computing with angular momentum and eccentricity vectors forces a 
redundancy as it resolves a degeneracy. In other words, in order to facilitate the traversal of high eccentricity and inclination 
regimes, we simulate a redundant 6A''-dimensional system of equations. The algorithm uses analytic averaging over the orbit of 
the perturbing particle (Gauss's method) and numerical averaging over the orbit of the perturbed particle. The timesteps are 
now fractions of secular time-scales, which are at least a factor of Mi,/{Nm) longer than steps taken in unaveraged simulations 
of equivalent accuracy. The price for this is that our method cannot capture the effects of mean-motion resonances, which 
appear in the secular dynamics at second order in the masses. 

The Gaussian ring algorithm provides us with a powerful numerical tool for the study of secular dynamics of nearly 
Keplerian N-body systems, especially those with large eccentricities and inclinations (so that perturbation expansions converge 
slowly or not at all) and/or crossing orbits. Our method resolves the dominant, secular evolution of these particles with much 
greater accuracy and speed than conventional N-body codes: in the examples we have examined the typical timesteps are two 
to three orders of magnitude larger, depending on the mass of the N-body cluster relative to the central mass (although the 
calculations per timestep are more lengthy). Moreover a single ring should be thought of as containing many stars, so the 
required for reliable results should be much smaller. 

We have discussed the implementation of the method O. The accuracy goals of the algorithm are encapsulated in two 
parameters, tquad and eint, the first expressing the desired fractional accuracy of the numerical quadrature of the forces over 
the perturbed ring, and the second expressing the desired accuracy of the integration of the equations of motion. The accuracy 
of the quadrature can be assessed by evaluating the rate of change of the semi-major axis ( jjS.lf) . which should be zero in 
secular dynamics. Provided one can live with moderate softening (b > O.Ola), between 16 and 128 force evaluations per ring 
should be sufficient to yield equad ~ 10~^^; 10 to 100 times that number of force evaluations may be required for extreme 
configurations or very small softening. 

Much can be done to enhance the efficiency and accuracy of the method: (i) adaptive quadrature routines can reduce 
the number of force evaluations required for ring pairs experiencing a close approach; (ii) geometric integrators for this 
most illustrious of Lie-Poisson systems should enhance the long-term accuracy of the numerical integration; (iii) the cost of 
quadratures makes the method eminently parallelizable in a distributed environment. 
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APPENDIX A: HALPHEN'S CONE: A DELAYED GEOMETRIC INTERLUDE 

As iHilll |l9o3) put it: 'Gauss first clearly indicated the role elliptic functions play in this subject. Halphen has since presented 
the investigation in an elegant manner.' Halphen's elegant approach is thoroughly geometrical, based on the realization that 
the average force exerted by an eccentric Keplerian ring at a given point depends only on t he invar i ants o f the cone formed 
by ring and point. Halphen's approach has received numerous renditions since; for example iMusenI (|l963l ) further simplifies 
Halphen's computations by noting the role played by a natural dyadic in the computation. Here we stick with the essentials 
of HiU's account, which is closest in spirit (and date) to Halphen (and Gauss). We describe the derivation for an unsoftened 
potential here, with the modifications required by softening in Appendix I All 

The origin of the reference frame is best taken at the attracted particle. Denote by r = {x, y, z) the attracting planet, 
and by r* = (a;*, t/*, z*) the central body, with r = \r\ (note the change in notation from the main text, where r denotes what 
is here r — r*). The averaged direct force is 

[A = |^/^, (Al) 



where I is the mean anomaly. Equal area in equal time implies d//(27r) = Aa / {-ko? \/l — e^), a being the area of a sector 
swept by the radius vector from the central body to the attracting ellipse, and 7ra^\/l — the area of the ellipse. Let 
w = (r — r*)xdr, where dr is the change in position of the attracting body over some small time interval > as it travels 
around its orbit. Then da = and h — iy-r*/|iu| is the perpendicular distance from the origin to the plane of the orbit 

of the attracting body. Thus da = ^iy-r*//i and 

, „,, Gm f r r^,-(r — r*) xdr Gm / rdr■r^,Xr , , „, 

= . , o <P - . , o n f <P -a . (A2) 



/C ^liiLU \ ± — ^ Jc 

where the curve C is the orbit of the perturber. 

Consider for example the a;-component of the force. Using Stokes's theorem equation l|A2|l may be re-written 

,r 1 Gm f , _ /xr*Xr\ Gm f , fixr^-r Xt\ , , „, 

\fx\i = , / ds-Vx — = , / ds-r ; , (A3) 
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where ds is an element of surface area and the integral is over the surface area bounded by the curve C, i.e., the area swept 
out by the vector r — in one orbit of the perturber. Now consider the cone having vertex at the origin (the attracted 
particle) that passes through the curve C, i.e., the set of all vectors parallel to r{l),0 ^ I < 2n. The curve C lies on this cone, 
and if it is smoothly distorted to some other curve on the cone the additional contribution to will be zero, because the 
normal to the cone ds at r is perpendicular to r. We conclude that the integral in equation (|A2p is independent of the curve 
C so long as it lies on the cone. We may therefore choose a curve C on the cone that simplifies the line integrals in (|A2|I . 

We write the equation of the cone in a reference frame with origin at the attracted particle and axes parallel to the x, y, 
z axes defined in ^2.1\ {x points from the central body to the periapsis of the orbit of the attracting particle; z is parallel to 
the angular momentum of this orbit, and y = zXx). In these axes the position of the central body is r* — Xi,x + j/*y + hz. 
The orbit of the attracting particle is described by the ellipse {1 — e^){x — ~\- ae)^ + {y — y*)'^ — {1 — e^) on the plane 
z = h. Thus the equation of the cone is 

(1 — e'^){xh — Xi,z + aez)^ + [yh — y-kz)^ — (1 — e^)a^z^ = 0, (A4) 
or 

r'^Wr = 0, (A5) 

where r — xx — yy + zz and 

/ (l-e^)/!^ (1 - e2)/i(ae - I*) \ 

W = ~yji . (A6) 

V (l-e^)/i(ae-a;*) -y^h {1 - e^){ae - x^^f + yl - [1 - e^)a^ ) 

This quadratic form can be diagonalized by an orthogonal transformation O to new Cartesian coordinates u — (iti, it2, Ms), 
i.e., r = Ou. Thus 

O'^WO = diag ifi,) (A7) 
and 

r'^Wr = diag {fii, ^2, lJ.3)u (A8) 

is the equation of the cone. The /^'s are the eigenvalues of W and the roots of the function 

= ^J.'^ + ~ e^)a'^ - yj - {2 - e'^)h'^ - {1 - e^){ae~ x^,f] 

+^Ih^{l - e^)[yj + h^ + (ae - x^f - (2 - e^)c?\ + (1 - e^fa^h". (A9) 

We note that 

2/(0) = (1 - e')'a'/i'' > 0, y{h^) = -ey.^h^ (AlO) 

which implies that the roots are real, one negative and two positive, and we may therefore assume that /xi > /X2 > > ^13. 
The eigenvectors of W may be written 

= [ 1) , (All) 

where 6k is real, chosen so that h^-h'' = 1. We assume that 61 and 62 are positive and then choose the sign of 63 so that ei, 
62, 63 form a right-handed triad, that is, 63 = ei Xe2. Then the orthogonal matrix O is defined by Oij = and detO = +1. 
The equation of the cone in the new coordinates is 

/iiui + ^2^2 = l^shs- (A-12) 
We choose the curve C to be the intersection of the plane u^, = const with the cone. A parametric equation for C is then 

ui = au-i cosijj, U2 = Pus sin 7/), a = v^j^sl/^ii, P = ^/\^J.3\/^J.2■ (A13) 
Then along C, 

= = Ui + U2 + Us = us{l + o? cos^ 1/) + sin^ i/") (A14) 
and 

dr-(rj, Xr) = du-(w* Xu) = Usd'>\)\Q.jiUi,z — a.u^,2 sin?/; — (iu^,\ cos 1/))]. (A15) 



The curve C as defined in equation (|A2|I is oriented counter-clockwise as viewed by a distant observer on the positive 
angular-momentum axis (the z axis). This is the same orientation as the curve defined by equation l|A13[) if and only if h 
and M3 (the distance of the integration plane from the origin) have the same sign. To account for this we assume 113 > and 
replace the pre- factor 27r/i in equation (|A2p by 27rj/i|. Then equation \K2\ becomes 
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, _ GmaP Z'^'" •tt3^J*3 - cos^ ^ - -1x20*2 sin^ tp ( ^^P^ 

^' ^ 2T:\h\a^VT^ Jo (l + a2cos2^ + /32sin2V)3/2 ' ^ ^ 

The integral can be evaluated using equations (|66|l . using the substitution -^/j = ^tt — T: 

[/'h = /(^^-^^ [(^Ga + G,)E{.) - (1 - .^)G.7^W1 (A17) 



Ga = 11*3U3 — lt*2U2, Gfc = lt*2W2 — lt*l Ul . (A18) 



where 

^2 ^ ImsI Ml - At2 
Ml M2 - M3 ' 

It is straightforward to demonstrate numericaUy that this expression for [/']; is the same as that given by equation (|67p . and 
that the argument k of the elhptic integrals in that equation is equal to their argument here, 



Al Halphen's cone, softened 

At first look, it appears that Halphen's geometric treatment, which is predicated on the homogeneity of the integrand, would 
fail when a length scale, the softening length, is introduced in the problem. However, the fact that Gauss's algebra is flexible 
enough to handle softened interactions suggests the opposite. So one wonders whether, despite softening, Halphen's derivation 
can be generalized to a softened potential. The answer is yes, and following are the essentials of how. 
With softening, equation (|A2|) is modified to 

where b is the usual softening length. We now stretch the coordinates in the directions x, y parallel to the orbit plane of the 
attracting particle, setting 

r — XX + yy + zz = sxx + syy + zz, r = xx + yy + zz. (A20) 
Then along any small increment of the curve C (the orbit of the attracting particle) 

dr-(r* Xr) = s^dr-(r* Xr). (A21) 
Furthermore on this curve 

r' + 6^ = + + /i^ + 6^ = s^{l£^ +y^) + + b\ (A22) 
If we now set 

5^ = 1 + 6^*^ then + fe^ = sV^ (A23) 
Thus equation (|A19|l becomes 

[/'], = [J^x+J^y + s-'Jj], where [/'], = \ (f ^dr-(r.Xr) ^ 



(A24) 



the contour C is obtained by shrinking the vectors from the attracted particle to the star and from the star to the attracting 
particle by a factor s in the directions parallel to the orbit plane of the attracting particle. 

Since the integrals for [/']; in equation HA19|I and for [/ ]; in equation (|A24[I are identical, the derivation in the preceding 
section can be carried through with only minor changes. The analog to equation (|A4|) for the cone is 

(1 — e^)(xh — Xi^zj s + aezj s)^ + (%jh — yi^zj sf' — (1 — e'')a' z^ j s' = 0. (A25) 
The cubic equation (|A9|I becomes 



= m"^ + M^[(l ~ e'^a' I s' — y*' I — (2 — e^)/i^ — (1 — e^)(ae — a;*)^/s^] 

^-Jih\\ ~ e')[y*Vs' + + (ae - x,f j - (2 - e^)a' j s'\ + (1 - e'fa'h^. (A26) 
The eigenvectors of W are 



(1 — e )h{ae~Xi,)/s y*h/s 



/ife - (1 - e2)/i2 ' h'^' 
The remaining formulae are the same, except that [/']; is to be replaced by [/'];, and fXi is replaced by Ji^ 



(A27) 
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APPENDIX B: GENERAL-RELATIVISTIC PRECESSION 

We may correct Newton's equations with the dominant secular contribution from general-relativistic effects. These effects are 
mainly feh by nearly radial orbits on their close encounters with the central black hole. When averaged over the orbital period 
these corrections contribute a Hamiltonian 

The gradients are V^iifgr — and 

^A^.--'-^^^- (B2) 
The secular equations of motion are l|Tremaine et al.ll2008l ) 

To this order, the main effect is on precession of periapsis, leaving eccentricity, orientation and semi-major axis of the orbit 
unaffected 

dA 67rGM* ^ . 

dF^ cW-e^P/^ '^^''- ^""'^ 

where P = 27ra^/^/(GM*)^''^ is the orbital period. 

For a spinni ng (Kerr) b lack hole, the dominant secular effect is Lense-Thirring precession. The relevant averaged Hamil- 
tonian is ()Landau fc Lifshitz 2005. ') 

^Ker. - -^^sj2(T^^S.L, (B5) 

where {GM^/c)S is the spin angular momentum of the central black hole and \S\ ^ 1. Then 

(,3^5/2 _ g2-j3/2 ^ g3a5/2(-X — g2)5/2 ^ ' ^ ' 



The equations of motion (|B3} yield 

^ = c3a3r2p(T-e2)3/2«X-^' ^ = ,3„3r2p(f 1^62)3/2 " 3I,(I,.5)] X TI. (B7) 



REFERENCES 



Fabrycky D., Tremaine S., 1997, ApJ, 669, 1298 

Ford E.B., Kozinsky B., Rasio F.A., 2000, ApJ, 535, 385 

Gauss C.F., 1866, Werke, 3, 331. Dieterich, Gottingen. See http://resolver.sub.uni-goettingen.de/purl7PPN235999628l (in 



Latin) 

Goriachev N.N., 1937, 'On the method of Halphen of the computation of secular perturbations'. University of Tomsk, 1-115 
(in Russian) 

Giirkan M.A., Hopman C, 2007, MNRAS, 379, 1083 

Hagihara Y., 1971, Celestial Mechanics II, Part I. MIT Press, Cambridge, MA., pp. 199-233 

Halphen G.H., 1888, Traite des Fonctions EUiptiques et de leurs Applications II. Gauthier-Villars, Paris 

HiU G.W., 1882, Astron. Papers Am. Eph. Naut. Aim., 1, 317 

HiU G.W., 1901, Am. J. Math., 23, 317 

Holman M., Touma J., Tremaine S., 1997, Nat, 386, 254 

Innanen K.A., Zheng J.Q., Mikkola S., Valtonen M.J., 1997, AJ, 113, 1915 

Jacobs v., Sellwood J., 2001, ApJ, 555, L25 

Kozai Y., 1962, AJ, 67, 591 

Landau L.D., Lifshitz E.M., 2005, The Classical Theory of Fields, 4th ed. Butterworth Heinemann, Amsterdam 

Laskar J., 1986, A&A, 157, 59 

Laskar J., Simon J.L., 1988, A&A, 157, 59 

Libert A.-S., Henrard, J., 2008, Celest. Mech. Dyn. Astron., 43, 37 
Mazeh T., Krymolowski Y., Rosenfeld G, 1997, ApJ, 477, L103 

Murray CD., Dermott S.F., 1999, Solar System Dynamics. Cambridge Univ. Press, Cambridge 



© 2008 RAS, MNRAS 000, 000-000 



28 J. Touma et al. 



Murray N., Holman, M., 1999, Sci, 283, 1877 
Musen P., 1963, Rev. Geophys., 1, 85 
Peiris H.V., Tremaine S., 2003, ApJ, 599, 237 
Rauch K.P., Tremaine S., 1996, New Astron., 1, 149 
Saha P., Tremaine S., 1992, AJ, 104, 1633 
Sambhus N., Sridhar S., 2000, ApJ, 542, 143 
Sridhar S., Touma J., 1999, MNRAS, 303, 483 
Touma J., 2002, MNRAS, 333, 583 
Touma J., Wisdom J., 1998, AJ, 115, 1653 
Tremaine S., 1995, AJ, 110, 628 

Tremaine S., 2001, Celest. Mech. Dyn. Astron., 79, 231 
Tremaine S., Touma J., Namouni F. 2008. larXiv:0809.0237l 
Wu Y., Murray N. 2003, ApJ, 589, 605 



© 2008 RAS, MNRAS 000, 000-000 



